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Abstract 

The aim of this study is to devise numerical methods for dealing with very high¬ 
dimensional Bermudan-style derivatives. For such problems, we quickly see that we can 
at best hope for price bounds, and we can only use a simulation approach. We use the 
approach of Barraquand &: Martineau [2] which proposes that the reward process should 
be treated as if it were Markovian, and then uses this to generate a stopping rule and 
hence a lower bound on the price. Using the dual approach introduced by Rogers [H] 
and Haugh & Kogan [8], this approximate Markov process leads us to hedging strategies, 
and upper bounds on the price. The methodology is generic, and is illustrated on eight 
examples of varying levels of difficulty. Run times are largely insensitive to dimension. 


1 Introduction. 

A general optimal stopping problem can be formulated as finding 

sup E[Zr] (1.1) 

0<t<T 

where the termination time T is a hxed positive real, the process Z is a given process adapted 
to some hltration (J^t)t>o, and r is constrained to be an (J^t)-stopping time. We shall in this 
paper assume that the reward process Z has the form 

= 9{ti Xt), ( 1 - 2 ) 

where X is some Markov process, and g is some measurable function of time and Xf. For 
what mainly concerns this paper - the analysis of Bermudan options - this formulation is 
already sufficiently general. Much of what follows is valid more widely, but at times we shall 
refer to properties that relate to a typical hnancial context; the aim is to provide bounds, and 
we need context in order to assess the effectiveness of these bounds. 

Of course, it is well known how to solve an optimal stopping problem of the form given 
by (HU), (HU; we dehne the value function 

V*{t,x)= sup E[g{T,X^)\Xt = x] (1.3) 

t<T<T 
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and find V* by dynamic programming. When this can be done, this is everything we could 
ask for. However, explicitly-soluble examples are rare, and we soon have to reach beyond 
them. For example, the standard Bermudan put option, where 

g(t,x) = — exp(x))~^ (1.4) 

and the process X is a Brownian motion, is a celebrated example where no closed-form solution 
is known. This has launched many studies in the last 25 years, and very efficient numerical 
schemes have been devised. But what if we have some vector = (X^, ..., Xf) of correlated 
Brownian motions, and the reward function is 

g(t,x) = e~^*(K — exp( minx^))~^ , (1.5) 

i 

a so-called min putl If we look for the value function, we have to determine numericalljH 
some function of d variables. The standard Bellman equation approach requires us to calculate 
recursively V)*(-) = H*(f, •) as 

V:{x) = max{ g{t,x), E[V,\,{X,+,)\Xt = x] } ; (1.6) 

but if d is large (fifty, say), how is to be stored? How is the expectation on the right- 
hand side of fll.bp to be calculated or approximated? The more one thinks about these issues, 
the clearer it becomes that calculating an approximation to the value function can only work 
in dimensions that are not too high, and will most likely rely heavily on the structure of the 
problem under study. In other words, any methodology that attempts to identify the value 
function will be of restricted applicability. 

So we must be content with less; but less may be enough. If we had calculated the value 
function, what use would we make of it? We would use it to determine the optimal control: 
at each time t, we would see the state Xt = x oi the system, and if Vf{x) > g{t, x) we would 
continue, otherwise we would stop. We would use it to determine a fair price to pay or ask 
for the derivative at time 0. We would use it to delta-hedge the derivative once it had been 
sold. The view taken in this paper is that the analysis of a Bermudan option requires just 
these: 


• At each time, whatever the state, we are able to decide whether or not to stop; 

• At each time, we are able to propose a hedge for the next period; 

• We are able to provide reasonably close bounds for the price of the derivative at all 
times. 

^For a Bermudan option, which has only finitely many possible times of exercise, the optimization is over 
a discrete set of times, though we generally think that the time of the underlying process runs continuously. 

^If we can only solve the one-asset problem numerically, we can certainly only solve the d-asset problem 
numerically. 
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Exact knowledge of the value function would achieve all of these objectives, but can we attain 
them without knowing the value function or an approximation to the value function? The 
message of this paper is that this can be done. We present an approach to Bermudan options 
with the properties: 

• The only information required about the Markov process is the ability to simulate a 
step of the process; 

• The methodology is generic - the same code is used to do the calculations for all exam¬ 
ples, changing only the specihcation of the Markov process X and the stopping reward 
9', 

• Computational cost is largely insensitive to dimension^ so a derivative written on hun¬ 
dreds of underlyings can be solved simply by changing the dimension parameter in the 
code, and takes only a little longer to run; 

• Upper and lower bounds on the price differ by typically 5-10 percent (sometimes more, 
sometimes less); 

• The stopping rules and hedging rules obtained are very simple to calculate and imple¬ 
ment. 

• Computational times depend on the problem, but a few tens of seconds usually suffices. 

The claim is that the method offered here is an effective general method for dealing with 
any Bermudan option. In truth, the components used are already known in one form or 
another, and what is added here is the judicious combination of them, and a redehnition of 
the questions to be answered. The waypoints are the following: 

• Any numerical scheme has to be a hnite calculation, so the Markov process has to be 
coerced to a hnite set of values; 

• The hnite coercion of the underlying Markov process X has to be tailored to the stopping 
reward - using the approach of Barraquand & Martineau [2] ; 

• The hnite coercion generates a stopping rule, and a hedging rule, using the dual approach 
of [H], [S], [I] which can then be evaluated by simulation. 

The particular structure of a problem may suggest variants that improve on the performance, 
but in the examples studied any improvement from simple variants is not large. 

2 The general methodology. 

The general situation concerns a Markov process X taking values in a statespace A, a stopping 
reward function g, and a hnite set T C [0,T] of size Nt of possible times to stoplf|. The aim 

^We will always assume that 0 and T are in T. 
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is to associate this problem with an optimal stopping problem for a (discrete-time) Markov 
process with a finite statespace, and the way this is done is simply a paraphrase of the method 
of Barraqnand & Martineau |2]. 

Definition 2.1. A (real-valued) Markovian coercion of X is specified by a measurable function 
: T X df I—)■ M, an Nt x A^bins matrix rj of bin values , an x (A^bins — 1) naatrix y of bin 
edges, and an {Nt — 1) x A^bins x A^bins array P with the properties: 


1. for each t G T, 




,(b 


.( 2 ) / .,( 2 ) 


vr < yfi' < vr < yr < ■ ■ < 


2. for each t G T\{T}, P{t, ■, •) is an A"bms x A^bins transition matrix. 

We can approximate the real-valued process g{t,Xt), t G T, by using the matrix y of bin 
edges to define a partition of M for each f G T: 

jb) = (_oo^ ^d)j ^ ^ jX) ^ qq) ('21) 

and then the matrix p of bin values to define the approximating process 

N 

k=l 

For the current purpose, we shall take g to be the function appearing in the dehnition fll.2p 
of the stopping reward. 

The process Y takes only finitely many values, but will not in general be Markov. Nev¬ 
ertheless, the essential idea of the Barraquand-Martineau approach is that we pretend that 
it is, with transition probabilities given by the array P. We then solve the optimal stopping 
problem for this Markov coercion, and use the solution found to propose exercise and hedging 
strategies for X. 


2.1 Implementation. 

For the implementation, we shall assume that the state process X takes values in some (subset 
of) euclidean space The dimension d of this space is unrestricted, and can be quite large. 
The numerical implementation consists of four stages: 

1. use simulation to initialize, calculating the transition matrix array P, and the bin edges 
and values; 

2. calculate the value and optimal stopping rule for the (presumed Markovian) process Y ; 
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3. find a lower bound by evaluating the performance of the stopping rule from step 2 when 
applied to stopping the process X; 

4. hud an upper bound by evaluating the hedging rule derived from step 2 when applied 
to the actual process X. 

We now give some more detail on each of these stages in turn. 

Initialization. First choose some number Wins of bins, and some number Wiock of points 
in each ‘half bin’, so that in total we will simulate Wim = 2WiockWins sample paths. The 
simulation will £11 up a Nt x d x Wim array, though since the process is Markovian, we can 
economize on space by just carrying along the current values, in effect a. dx Wim array, which 
gets overwritten with newly-calculated values. 

So suppose that we have the d x Wim array x of the valued at time ti G T. We now 
use the simulation of the Markov process to step these values forward to the next time fj+i, 
creating a d x Wim array x'. Next apply the function g to give 

yj = g{ti,x[:J]), Vj = g{ti+i,x'[:,j]) (j = 1,..., Wim)- 

If we define ?/(i) < j/( 2 ) < ... to be the sample (yj) in increasing order, we define the bin 

edges at time ti to be the values y(2kNuoA)) ^ = 1) • • • > Wins — 1, and the bin values at time 
ti to be the values l/((2fc-i)Vbiock); k = 1,..., Wins- We similarly calculat^ the bin values and 
edges at time fj+i. Now for each simulated path j we can see which bin that path was in at 
time ti, and which bin that path moved into at time tj+i; counting the number of paths which 
moved from bin £ at time ti into bin m at time tj+i gives us an estimate of the transition 
probability P[i,£,m]. 

It is worth remarking that this step differs slightly from what Barraquand & Martineau 
do; we let the data tell us where the bin edges should be, and Barraquand & Martineau set 
the bin edges before any simulation takes place, as an a priori modelling choice. A similar 
use of the simulated data to determine an approximation procedure appears in Bouchard & 
Warin |1]. 

Calculating the value function and optimal stopping rule. Now that we have made 
a Markov chain proxy which jumps at discrete times from one bin to another according to 
the transition probabilities stored in P, the calculation of the value V and stopping rule S is 
done by dynamic programming. The arrays V and S are both x Wins ; the values in S 
are either 0 or 1, where S[i,k] = 1 signifies that we should stop if at time ti the Wvalue is 

'*In accordance with Python notation, we denote the {i,j) element of an array z by z[i,j], and the jth 
column of z by z[:,j]. 

®We implicitly assume that the y-values are distinct. This is clearly not going to happen for (say) a put 
option, but we make this happen by replacing the reward maxjO, K — S} by the reward maxfefiy —S'), (K — S)} 
for some small e. The error committed will be small compared to other errors. Likewise, we do not bother to 
locate the bin edges in between values of y(j) as perhaps we ought. 
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in bin k. The time taken to compute V and S is negligible. As a notational convenience, for 
?/ G M we shall write 

y{t,y) = v[i,k], S{t,y) =S[i,k] (2.3) 

(k) 

when t = ti eT and y E , the fcth bin at time ti. 

Lower bound. The computed stopping rule S is now used to provide a lower bound for the 
value of the option. We simply simulate a large number of paths of the process X, and for 
each path we stop the hrst time U for which S(ti, g{ti, Xt^)) = 1. For each path, this gives a 
stopping value, and the lower bound is just the average over all paths. 

Upper bound. To derive an upper bound, we need to recall some results about the dual 
approach from 111 . 0 . where it is shown that the value fll.ip of the optimal stopping problem 
has the alternative characterization 


sup E[Zr]= min E[ sup (Zt — Mt)], (2.4) 

o<T<r MeMo o<t<T 


as a minimum over the space Ado of martingales vanishing at 0. The interpretation of any 
M G Ado as a hedging martingale is explained in HI; this is an important part of the 
approach adopted here, because the approximate hedge is achieved by a ‘good’ martingale 
from the dual approach, not by delta-hedging. As a consequence, it is not essential that we 
can evaluate an option price to high accuracy; the main driver for demanding high accuracy in 
pricing is to perform delta hedging, yet the present approach completely steps around all such 
considerations. Indeed, the entire method presented here works for any Markov process X, 
including finite-state Markov chains, for which the concept of delta hedging is as meaningless 
as differentiating with respect to an integer argument. 

The minimum on the right-hand side is attained when M is the martingale part of the 
Snell envelope process: see [2] . In this setting, the Snell envelope process is simply the value 
function V*{t,Xt) evaluated along the path, so the optimal martingale difference sequence 
would be 


- K = I j-j. 


(2.5) 


Of course, we do not know V*, but we have calculated and stored in the array V some approx¬ 
imation to U*, so a natural approximation to M* would be found by taking the martingale 
difference sequence 


M,,„ - M,, = V(t,+u Y,,„) - E\V(U+u y;,„) I T,,], (2.6) 

where Y is as dehned at fl2.2p . This is essentially the approach of Andersen & Broadie [1]. The 
only issue with this is how we are to evaluate the conditional expectation on the right-hand 
side of 02.61) . The function U(ti+i, •) is a simple function, taking only Xbms values, but if we 
have simulated a sample path and we see X^, = x, how are we to calculate (or approximate) 
£'[U(tj+i, Ui+i) I Xt. = X \I The approach adopted is to perform a subsimulation of some 
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Nsuh values of starting from Xt. = x. It is generally considered a bad idea to perform 
subsimulations, because this will usually take a lot of time, and may not be very accurate, 
but in this application the approach is effective because is scalar. The importance of this 
is that our subsimulations do not need to search out some high-dimensional space, they only 
need to search out the real line; in practice, if the payoff g is continuous, and the time-step 
not too large, most values of will be fairly close to gift, XtJ, so a relatively small 

number of subsimulations will suffice. In the examples reported later, we used only a few tens 
of subsimulations, usually on 4000 simulated paths; full details are reported for each example 
studied. 

This explains how we construct candidate hedging martingales from the array V. In a 
hnancial context, we would want to be able to express such martingales in terms of traded 
assets; this could be done by calculating delta-hedges for the martingale differences arising 
from the approximation, but since this would be application-specihc, we prefer not to go into 
the detail of how this would be done. The approach offered merely indicates a martingale to 
be used for hedging, not how exactly this is to be synthesised from marketed assets. 

3 Examples. 

Here we present a range of examples to illustrate the methodology, some familiar from the 
literature, others not. In the following examples, we will be considering a d-vector of log- 
Brownian assets whose prices SI at time t evolve as 

dSi = Si (3.1) 

V j=i / 

where W is a. d-dimensional Brownian motion, the aij and /ij are previsible processes which 
are assumed constant, except in Example 13.81 We write 

$ =aa^, (3.2) 

a positive-dehnite symmetric matrix. Since the focus is on derivative pricing, we shall assume 
that we are working in the risk-neutral measure, which amounts to the condition 

g.i = r-\%ii, (3.3) 

where r is the riskless rate of interest, assumed constant, except in Example fl3.8p . 

It proves convenient in most of the examples to simulate the discounted log prices 

xl = - f rsds + log SI = ^ (TijW^ - i Sjj t log S'q. (3.4) 

Jo ^ 
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Variants. If V is a positive martingale, Nq = 1, then 


snp E 

0<t<T 


9{r,Xr) 


sup E 

0<t<T 



9{r,Xr) 

Nr 


sup E 

0<t<T 



9{r,Xr) 

Nr 


sup E 

0<t<T 


9{'r,Xr) 

Nr 




where 

§ = Nr ( 3 . 5 ) 

If the martingale N can be expressed as N^ = X^), then we can work in the new measure 

P with the new reward g{t,x) = g{t,x)/'ilj{t,x), and it may be convenient sometimes to do 
this. This change of numeraire approach is reminiscent of Jamshidian’s [9] version of dual 
American option valuation. We could similarly transform the stopping reward g{t,Xt) to 
g{t, Xt) — (fit, Xt) + (p(0, Xq), where (p is some function for which (p(f, Xt) is a martingale, for 
example, the value of the European option. 


In the examples which follow, we calculate the simulation-based upper and lower bounds 
for various derivative prices. In the tables presented, we report also the simulation values 
of the European equivalent derivative. This should of course always be cheaper than the 
Bermudan option. In some places this natural inequality appears to be violated by a small 
amount. Partly this is because the standard error inherent in the simulation method (reported 
in parentheses after each entry in the tables), but also because the European prices are 
computed from the averages of the actual values of the option at expiry, whereas the lower 
bounds are based on the discretized values of the the option at expiry, and these are not the 
same. Discrepancies are in any case all very small even when present. 

The total compute time is reported also; usually, the main part of the calculation was 
the calculation of the dual upper bound. Sometimes the initial calculation of the transition 
matrices was quite time consuming, but typically less than the upper bound. 


3.1 Min put. 

In this example, the state variable is Xt = Xt, the vector of discounted log prices, and the 
reward function for stopping at time t will be 

g{t,Xt) = ( Xe“''* - exp( min X;) )^. (3.6) 

l<2<a 

This example was studied in ng, and the hgures in the column MC price of Tabled] were taken 
from that paper; for d = 30, 60 no values are given in H. The column headed European is 













the simulation value for the European option, where no early exercise is allowed. The columns 
headed low and high are sample means obtained from a simulation method. Standard errors 
are reported in brackets after the mean values. Notice that the lower bound is in all cases less 
good than the bound from the European price. This may be in part due to the hnite sizes 
of the bins, and the error arising from that, but it indicates that this simple approach is not 
able to extract the holder’s early exercise value of the option. Nevertheless, the bounds are 
reasonably close for all values of d, even quite good for larger values of d, and the run times 
are going up roughly linearly with d. 


d 

European 

low 

MC price 

high 

gap(%) 

time 

2 

24.78 (0.07) 

24.71 (0.08) 

25.16 

25.65 (0.30) 

3.40 

8.89 

3 

31.27 (0.06) 

31.16 (0.08) 

31.76 

32.29 (0.34) 

3.16 

10.20 

4 

35.75 (0.06) 

35.74 (0.07) 

36.28 

36.81 (0.35) 

2.88 

11.42 

5 

39.22 (0.06) 

39.12 (0.07) 

39.47 

39.84 (0.33) 

1.56 

12.49 

10 

48.01 (0.04) 

47.81 (0.05) 

48.33 

48.50 (0.28) 

1.00 

18.64 

15 

52.13 (0.04) 

52.06 (0.04) 

52.14 

52.62 (0.23) 

0.93 

26.38 

30 

57.82 (0.03) 

57.66 (0.03) 

- 

58.21 (0.17) 

0.66 

61.20 

60 

62.28 (0.02) 

62.18 (0.03) 

- 

62.57 (0.16) 

0.47 

183.49 


Table 1: Min put prices. The d assets are independent, S'*(0) = 100. Other parameters 
are K = 100, r = 0.06, T = 0.5, an = 0.6. Parameters for the simulations are N"bms = 
200, iVbiock = 200, Nt = 40, iVprimai = 50000, Nduai = 400, Nsub = 60. 


3.2 Max call. 

As with the min put, the state variable is Xt = Xt, the vector of discounted log prices, but 
this time the reward function for stopping at time t will be 

5f(t,XJ = ( exp(m^Xj) (3.7) 

In order to make the problem interesting, the assets will be assumed to pay dividends at a 
constant rate. This example was studied by Broadie & Glasserman^ [6], and has been used as 
a test example in a number of other studies, including 0 . 0 , 0 . [To] . Notice that this time 
the lower bound we obtain is signihcantly bigger than the European price, so the holder is 
able to use the simple Markov coercion heuristic to extract some of the early exercise value. 

®In the preprint version, the example presented in Table [2] is said to give results for common expiry T = 3 
with 3, 6, 9 equally-spaced exercise opportunities, but the published version gives the same numerical values 
purportedly for end-of-year exercise opportunities with expiries 3, 6, 9 years. Our numerics show that in fact 
the problem solved is correctly stated in the preprint, and mis-stated in the published version. 
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The bounds are less close than for the min put example, but this is a more difficult option to 
handle; errors of the order of 5-6% need to be improved, but already give usable information. 


m 

^0 

European 

low 

BG price 

high 

gap(%) 

time 


90 

14.62 (0.06) 

15.39 (0.08) 

16.006 

16.18 (0.10) 

4.87 

4.40 

3 

100 

22.98 (0.08) 

24.37 (0.09) 

25.284 

25.43 (0.12) 

4.15 

4.38 


110 

32.60 (0.09) 

33.93 (0.11) 

35.695 

35.86 (0.14) 

5.39 

4.49 


90 

14.58 (0.06) 

15.92 (0.07) 

16.474 

16.65 (0.08) 

4.39 

8.00 

6 

100 

23.14 (0.08) 

24.95 (0.08) 

25.290 

26.31 (0.10) 

5.17 

8.35 


no 

32.64 (0.09) 

35.03 (0.09) 

36.479 

37.24 (0.13) 

5.94 

8.30 


90 

14.57 (0.06) 

16.05 (0.07) 

16.659 

16.93 (0.09) 

5.20 

11.77 

9 

100 

23.05 (0.08) 

25.14 (0.08) 

26.158 

26.88 (0.11) 

6.47 

12.08 


no 

32.61 (0.09) 

35.23 (0.09) 

36.782 

37.57 (0.12) 

6.23 

11.83 


Table 2: Max call prices on 5 independent assets with common volatility a = 0.2 and expiry 
T = 3. There are m = 3,6,9 exercise opportunities at times zT/m, i = 0,... ,m. Other 
parameters are K = 100, r = 0.05, 5 = 0.1. Simulation parameters are iVbins = 500, A^biock = 
100, iVprimal = 50000, A^dual = 4000, N^nh = 150 . 


3.3 Basket put. 

This is an example studied in Kovalov, Linetsky & Marcozzi m, and subsequently in Jin et 
al [To]. The state variable is again the vector of d discounted log prices, and this time the 
stopping reward function is 

d 

9(i, X,) = ( A'e-" - d-^ (3.8) 

i=l 

All the stocks start at 100, the strike is K = 100, the riskless rate is 0.03, the expiry is 
T = 0.25, and the individual asset volatilities are all 0.2, but this time the assets are not 
supposed independent; there is constant correlation p = 0.5 between all the asset^. Kovalov 
et al use a numerical PDE approach, Jin et al use a simulation methodology, and both 
approaches appear to give better precision than the method we have used here. Nevertheless, 
as we shall soon see, the difference in precision is not practically relevant. We report the 
results in Table [21 

These values were computed assuming that the volatility parameter a is equal to 20%. 
But are we sure of that? In any application, the volatility (assumed constant) would have to 

^The specification of the various parameters in m is internally inconsistent, and not in agreement with 
the parameters quoted in m- For the reported prices to be correct, we find that the parameter values used 
must be those we have stated. 
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d 

European 

low 

KLM 

high 

gap(%) 

time 

2 

3.08 (0.01) 

3.13 (0.02) 

3.14 

3.25 (0.01) 

3.59 

41.89 

3 

2.89 (0.01) 

2.93 (0.02) 

2.94 

3.04 (0.01) 

3.64 

48.03 

4 

2.78 (0.01) 

2.81 (0.02) 

2.84 

2.94 (0.01) 

4.35 

52.39 

5 

2.72 (0.01) 

2.75 (0.02) 

2.77 

2.87 (0.01) 

3.99 

57.54 

6 

2.68 (0.01) 

2.73 (0.01) 

2.73 

2.83 (0.01) 

3.81 

63.42 

12 

2.56 (0.01) 

2.60 (0.01) 

- 

2.70 (0.01) 

3.56 

99.03 


Table 3: Basket put. All assets start at 100, K = 100, T = 0.25, r = 0.03. All assets 
have volatility 20%, and the correlation between assets is 0.5. Other parameters are Abins = 
500, Abiock = 200, Nt = 40, Aprimai = 50000, Aduai = 1000, Asub = 160. 


be estimated; are we really sure that the volatility is not 19%? Or 21%? Are we really sure 
that the volatility will remain constant at 20% until expiry of the option? Suppose we repeat 
the calculations of Table [3] for those values of the volatility parameter and see what ranges 
for the price result. The outcomes are recorded in Table IH and what we see is that there is 
no overlap between the computed intervals for the price for the three values of a. In other 
words, the uncertainty in the price arising from our simulation bounds is comparable to the 
uncertainty in price which would arise from the uncertainty in the input parameter values. 


d 

low (19%) 

high (19%) 

low (20%) 

high (20%) 

low (21%) 

high (21%) 

2 

2.98 

3.07 

3.13 

3.25 

3.28 

3.42 

3 

2.76 

2.87 

2.93 

3.04 

3.09 

3.20 

4 

2.68 

2.79 

2.81 

2.94 

2.97 

3.10 

5 

2.61 

2.72 

2.75 

2.87 

2.93 

3.04 

6 

2.58 

2.68 

2.73 

2.83 

2.86 

2.97 


Table 4: Basket put. Parameters as for Table [21 except the volatility parameter which takes 
values 19%, 20%, 21%. 


3.4 Fixed strike Bermudan-Asian call. 

In this example, there is a single asset S, and the reward for stopping at time r is 

g{T,X,)=e-^-{A,-Ky, (3.9) 


where we define the average price 


A+ — 


Su du 
t + 5 


(3.10) 
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Here, 5 > 0 is some initial window required to prevent wild oscillations. There is also some 
initial lock-out time t* > 0 during which exercise of the option is forbidden. The state variable 
of the problem is Xt = [St,At,t]. The numerical results are reported in Tabled The gaps 
between the bounds are quite variable. In fact we used three different numeraires (the bank 
account, the stock, and the martingale Et[AT])- The results are usable, but not particularly 
good, and this is really because this example is intrinsically two-dimensional, so any attempt 
to coerce it to one dimension is missing something essential. We cannot expect a stopping rule 
which only looks at At to do very well, because the current value of St has to be considered 
as well; if St is high enough, the value g(t,Xt) is actually increasing, so we would certainly 
not stop at such a time. But a rule that only considers At would not understand that. 


^0 

^0 

European 

low 

FD price 

high 

gap(%) 

time 


80 

0.956 (0.008) 

0.945 (0.012) 

0.949 

0.974 (0.033) 

1.83 

29.04 


90 

3.230 (0.016) 

3.216 (0.023) 

3.267 

3.374 (0.073) 

4.27 

29.51 

90 

100 

7.571 (0.025) 

7.568 (0.035) 

7.889 

8.048 (0.116) 

5.93 

29.39 


no 

13.78 (0.03) 

13.77 (0.04) 

14.538 

14.81 (0.16) 

6.97 

29.72 


120 

21.21 (0.03) 

21.00 (0.05) 

22.423 

22.94 (0.18) 

7.58 

30.01 


80 

1.098 (0.009) 

1.088 (0.013) 

1.108 

1.160 (0.036) 

5.36 

29.39 


90 

3.557 (0.017) 

3.578 (0.024) 

3.710 

3.759 (0.080) 

4.82 

29.40 

100 

100 

8.133 (0.025) 

8.148 (0.036) 

8.658 

8.905 (0.127) 

8.50 

29.40 


no 

14.73 (0.04) 

13.92 (0.05) 

15.717 

16.39 (0.17) 

10.12 

29.62 


120 

22.09 (0.05) 

21.88 (0.07) 

23.811 

24.46 (0.21) 

9.67 

27.32 


80 

1.260 (0.011) 

1.233 (0.015) 

1.288 

1.337 (0.048) 

5.75 

30.32 


90 

3.911 (0.023) 

3.978 (0.033) 

4.136 

4.506 (0.107) 

11.72 

27.19 

no 

100 

8.885 (0.029) 

8.378 (0.040) 

9.821 

10.822 (0.133) 

17.90 

29.64 


no 

15.53 (0.04) 

15.61 (0.05) 

17.399 

18.40 (0.16) 

15.14 

29.69 


120 

23.02 (0.05) 

23.36 (0.07) 

25.453 

26.42 (0.20) 

10.92 

26.94 


Table 5: Fixed strike Bermudan Asian call. Parameters are a = 0.2, K = 100, t* = 0.25, 
5 = 0.25, T = 2. Other parameters are Wins = 500, Wiock = 100, Nt = 40, Aprimai = 
50000, Wuai = 4000, Wub = 125. 


3.5 Floating strike Bermudan Asian call 

The story is very similar to Section 13.41 except that the reward is 

g{T,Xr)=e-^^{Ar-Sr)^ (3.11) 

for stopping at time r. This example is in fact much easier than the hxed strike, because the 
process At/St is a Markov process already. Scaling says we need only vary 5*0 while keeping Aq 
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fixed, which is what we do. The results obtained when we £x Aq = 100 are given in Table |H1 
Once again, the gaps between the two bounds are usably close, and the lower bounds are well 
clear of the European values, so here the approximate stopping rule gives a very substantial 
improvement. 


^0 

European 

low 

high 

gap(%) 

time 

80 

3.98 (0.04) 

10.62 (0.03) 

10.81 (0.07) 

1.77 

8.88 

90 

3.936 (0.041) 

8.347 (0.032) 

8.618 (0.051) 

3.14 

9.21 

100 

3.866 (0.042) 

7.136 (0.033) 

7.485 (0.040) 

4.66 

9.11 

110 

3.873 (0.043) 

6.558 (0.035) 

6.909 (0.043) 

5.09 

8.93 

120 

3.942 (0.046) 

6.137 (0.037) 

6.529 (0.036) 

6.00 

8.82 


Table 6: Floating strike Bermudan Asian call. Parameters are a = 0.2, Aq = 100, t* = 0.25, 
6 = 0.25, T = 2. The calculations were done in the numeraire of the discounted asset 
price. Other parameters are A^bins = 200, A^biock = 100, Nt = 40, Aprimai = 50000, A^duai = 
4000, A'sub = 50 . 


3.6 Fixed window lookback option. 

This example illustrates the capacity of the methodology to handle high-dimensional prob¬ 
lems. Here we suppose that the stock price is recorded at times which are multiples of some 
h > 0, and stopping at time t = kh delivers reward 

g{T,X^)= sup Sjh, (3.12) 

k—a<j<k 

where a is some positive integer. This time, the state variable X has to record the last a 
values of the price, since the sup and inf are taken over a hxed window. The results are 
presented in Table [71 Notice that for this example the bounds are very close, getting slightly 
less good as the lookback parameter a rises. Equally noteworthy is the fact that the run times 
are changing little as we increase the lookback parameter; so in the hnal row of the table, the 
state variable is 25-dimensional. It should not be a surprise that there is so little variation in 
run times, because increasing a makes no difference to the simulation load; each period, we 
simulate one new value, all that is different is that we are storing more or fewer values from 
the past. 

3.7 Fixed window range option. 

This example is similar to the fixed window lookback option of Section 13.61 except that the 
reward for stopping at time t = kh is 

g{T,Xr)= sup Sjh- inf Sjh, (3.13) 

k—a<j<k k CL'^j'^k 
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a (days) 

European 

low 

high 

gap(%) 

time 

5 

103.50 (0.02) 

117.67 (0.01) 

118.51 (0.03) 

0.71 

93.00 

10 

105.96 (0.03) 

124.71 (0.02) 

126.46 (0.05) 

1.38 

88.65 

15 

107.80 (0.04) 

129.24 (0.03) 

131.86 (0.06) 

1.98 

91.34 

20 

109.38 (0.05) 

132.62 (0.03) 

135.74 (0.07) 

2.30 

96.68 

25 

110.79 (0.06) 

135.35 (0.04) 

138.77 (0.07) 

2.47 

100.58 


Table 7: Fixed window lookback option. Parameters were T = 1, a = 0.5, r = 0.05, 
Sq = 100, and the time interval was divided into 500 equal (half-day) time steps. The 
calculations were done in the numeraire of the discounted asset price. Other parameters were 
Tobins = 250, iVblock = 60, Nt = 500, lYprimal = 50000, A^dual = 4000, A^"sub = 25. 


where a is some positive integer. The results are reported in Table [HI The gaps between the 
upper and lower bounds are higher than for the lookback example expressed as a percentage, 
but the arithmetic gaps are roughly comparable, with similar run times. 


a (days) 

European 

low 

high 

gap(%) 

time 

5 

9.94 (0.03) 

23.97 (0.02) 

25.35 (0.04) 

5.43 

31.56 

10 

16.77 (0.05) 

33.31 (0.03) 

36.63 (0.07) 

9.05 

35.85 

15 

21.94 (0.06) 

39.56 (0.03) 

44.63 (0.09) 

11.36 

38.05 

20 

26.25 (0.07) 

44.36 (0.04) 

51.00 (0.11) 

13.01 

42.28 

25 

30.04 (0.08) 

48.39 (0.04) 

56.34 (0.14) 

14.12 

45.00 


Table 8: Fixed window range option. Parameters were T = 1, a = 0.5, r = 0.05, Sq = 100, 
and the time interval was divided into 250 equal time steps. The calculations were done in 
the numeraire of the discounted asset price. Other parameters were A^bins = 200, A^biock = 
50, Nt = 250, = 50000, iVa^ai = 4000, iV^ab = 25. 


3.8 Min puts with stochastic volatility and interest. 

In this example, we consider a situation where there are d > 1 assets, with stochastic volatility 
and interest rates. There are examples of such models in various places in the literature, for 
example, Medvedev &: Scaillet [13], Boyarchenko & Levendorskii |5], Jin et al. [10]. Heston 
dynamics for the asset and the volatility are popular in theoretical work, but it is far from 
clear that the dependence of the volatility of volatility on level takes the square-root form 
postulated in the Heston model; and still less is it persuasive that the Cox-Ingersoll-Ross 
interest-rate model correctly describes the volatility of interest rates when those rates are 
low. As this is a simulation study, we are freed from any need to choose models that are 
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theoretically tractabl^, so we may make some modelling assumptions that match observed 
behaviour better. 

So our story supposes that there is some market Brownian motion , and that log prices 
xl = log SI evolve as 

dxl = <T,‘{ PS dwf + p's dwf'' } + (r, - 1 (a<,)-‘) dt (3.14) 

where r is the riskless rate process, and ps € (0,1) is the correlatiorHof the log prices with the 
market Brownian motion. The process is a d-dimensional Brownian motion independent 
of W^. 

The volatility process a appearing in fl3.14p is represented as 

a* = d*exp(^*) (3.15) 

in terms of constants d* and a process ^ which is an OU process evolving as 

dt + (pg dWl^ + p' dWf ). (3.16) 

Finally, our model for the interest-rate process r is just a Black-Karasinski model: we 
have Tt = fexp{zt), where 


dzt = -I3r dt + ar {pr dW^ + p; dWl) (3.17) 

for some constants > 0, > 0 and p^, which would typically be assumed positive since 

we expect that as the market rises the rate of interest should also rise. 

Altogether then, this is a simple but sprawling model; even assuming (as we do here) that 
correlations are common across stocks, the parameter vector is 

d = {ps,Pi,Pr, (d*),f,/3^,cr5,/3^,a^). (3.18) 

What would be reasonable values for these parameters? For the interest rate evolution, we 
shall be guided by Black & Karasinski [3] and take f = 0.06, cXj, = 0.12, (3^ = 0.02, and 
Pr = 0.3. Correlations between stocks are variable, but typically in the range 0.25-0.60; we 
shall take ps = 0.3. For simplicity we assume all stocks have common volatility d* = 0.6. 
Fluctuations in volatilities are of the order of tens of percent, so by comparing with the 
standard deviation of an OU process, we impose 

^ = 0.1. (3.19) 

®Tractability is in any case illusory; we regard a model as tractable if there is a closed-form solution for a 
small number of derivative prices, overlooking the fact that for the majority of derivative prices there is no 
closed form solution. 

®When p is a correlation coefficient, we use p' to denote yl — fp. 
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We shall set (3^ = 4.5, so that the mean reversion time for volatility is of the order of three 
months, and this gives from fl3.19p that = 0.3. Finally, we take = 0.3. 

For the rest of our discussion, we shall focus on the case where d = 5, that is: there are 
hve assets. We will also restrict attention to min puts, taking strike K = 100 and all assets 
starting at 100 throughout. We shall also assume that the expiry T = 0.5. This will allow us 
to investigate the effects of varying initial values of r and cr, as well as various parameters. If 
we make ps = 0 and make the mean reversion parameters f3r, (3^ very large, we have in effect 
got back to the situation with independent assets, constant volatilities and interest rate that 
was studied in Examples 13.1113.21 So we should see the same answer; and we do - the range 
from this calculation came out to be [39.06, 40.65], which is the same as we found in Tabled! 
to within sampling error. 

Next, we can see what happens when we keep the volatility and interest rate constant, 
but allow ps to vary; the results are in Table O What we see is that while the correlation is 
not too far from zero, there is no clear effect on the price, but as the correlation between the 
assets rises, the price of the min put falls. This is to be expected; the higher the correlation 
the more alike the assets are, so there will be less dispersion in the prices at any time, so the 
minimum will be higher. 


Ps 

European 

low 

high 

gap (%) 

time 

-0.15 

38.81 (0.11) 

38.81 (0.07) 

39.95 (0.16) 

2.85 

19.87 

0.00 

39.32 (0.11) 

39.09 (0.07) 

40.16 (0.16) 

2.09 

20.10 

0.15 

38.85 (0.11) 

38.72 (0.07) 

39.90 (0.16) 

2.63 

20.11 

0.30 

37.58 (0.12) 

37.63 (0.07) 

38.87 (0.17) 

3.17 

20.16 

0.45 

36.00 (0.13) 

35.82 (0.08) 

37.06 (0.19) 

2.88 

20.15 

0.60 

33.33 (0.13) 

33.00 (0.08) 

34.15 (0.20) 

2.40 

20.12 


Table 9: Prices of min puts as ps varies. Volatility is constant at a = 0.6, interest is constant 
at r = 0.06. Other parameters are Wins = 200, Wiock = 50, W = 40, Wrimai = 50000, Wuai = 
4000, Wub = 50. 


Now we relax the assumption of constant interest rate, and let the interest rate evolve as 
in the Black-Karasinski specihcation fl3.17p . hxing ps at its default value 0.3, and observing 
the effects of different initial values of the riskless rate. We hold the volatility constant at the 
default value 0.6. The results are reported in Table fTOl As the initial interest rate risei^. the 
price of the min put falls, as would be expected from the stronger discounting of the stopping 
reward; even allowing for the fact that we can only obtain an interval for the price, the effect 
of change of initial interest rate is discernible. 

Having seen the effect of changing initial interest rate while volatility is held constant, we 
next hold the interest rate constant at its default value f = 0.06, and allow the volatility to 

^°The row of the table for rp = 0 was obtained by taking logrp = —15.6. 
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ro 

European 

low 

high 

gap (%) 

time 

0.00 

40.61 (0.12) 

40.45 (0.07) 

41.54 (0.18) 

2.26 

19.64 

0.025 

39.34 (0.12) 

39.27 (0.07) 

40.16 (0.18) 

2.03 

19.58 

0.06 

37.76 (0.12) 

37.50 (0.07) 

38.70 (0.17) 

2.42 

20.04 

0.10 

35.88 (0.12) 

35.86 (0.07) 

36.99 (0.17) 

2.99 

19.68 


Table 10: Prices of min puts as tq varies. Volatility is constant at a = 0.6, correlation between 
assets is ps = 0.3. Parameters are pr = 0.3, f = 0.06, f3r = 0.02, cxr = 0.12. Other parameters 
are Vbms = 200, iVbiock = 50, Nt = 40, Vpnmai = 50000, iVduai = 4000, Ngnh = 50. 


be stochastic. The impact of different initial levels of volatility is shown in Table [HI where 
it is assumed that the initial volatility is common across all the assets. Again, we see quite 
pronounced effect of the initial volatility on the price of the min put option; the price rises as 
the initial volatility increases, as one would expect. 


o-Q 

European 

low 

high 

gap (%) 

time 

0.10 

22.17 (0.09) 

22.14 (0.05) 

22.89 (0.09) 

3.15 

19.67 

0.20 

26.78 (0.10) 

26.66 (0.06) 

27.42 (0.12) 

2.34 

19.98 

0.30 

30.08 (0.10) 

29.99 (0.06) 

31.03 (0.13) 

3.04 

19.97 

0.40 

32.83 (0.11) 

32.86 (0.07) 

33.74 (0.15) 

2.59 

20.08 

0.50 

35.37 (0.11) 

35.33 (0.07) 

36.40 (0.16) 

2.84 

20.09 

0.60 

37.49 (0.12) 

37.50 (0.07) 

38.65 (0.17) 

2.97 

20.10 


Table 11: Prices of min puts as cxo varies. Interest is constant at f = 0.06, correlation between 
assets is ps = 0.3. Parameters are p^ = 0.3, d* = 0.6 for all i, (3^ = 4.5, = 0.3. Other 

parameters are Abins = 200, Abiock = 50, Nt = 40, Apnmai = 50000, Aduai = 4000, N^nh = 50. 


The hnal study considers the full model where both volatility and interest rate are stochas¬ 
tic. There are too many parameters to explore in a paper, so we content ourselves with holding 
the parameters hxed at their default values, and varying the initial riskless rate and initial 
volatility. We see the results in Table [Hi Again the comparative statics behave as one would 
expect, with the magnitudes of the effects of changing initial values being big enough to show 
up clearly, even allowing for the fact that we have only got bounds on the prices. 

4 Conclusions and discussion. 

The aim of this paper has been to see to what extent we are able to solve the problem of 
pricing Bermudan options in very high dimensions. Once we accept that for such problems it 
is impossible that we can know the value function, we realize that in fact various approaches 
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O-Q 

European 

low 

high 

gap (%) 

time 


0.2 

29.71 (0.10) 

29.53 (0.06) 

30.38 (0.12) 

2.18 

19.47 

0.00 

0.4 

35.77 (0.11) 

35.65 (0.07) 

36.45 (0.15) 

1.86 

19.64 


0.6 

40.38 (0.12) 

40.36 (0.07) 

41.34 (0.17) 

2.32 

19.73 


0.2 

26.85 (0.10) 

26.69 (0.06) 

27.61 (0.12) 

2.74 

19.84 

0.06 

0.4 

33.00 (0.11) 

32.78 (0.07) 

33.60 (0.15) 

1.79 

20.20 


0.6 

37.67 (0.12) 

37.70 (0.07) 

38.66 (0.17) 

2.49 

20.00 


Table 12: Prices of min puts as ao and tq vary. Parameters are Ps = Pi = Pr = 0.3, 
d* = 0.6 for all i, f = 0.06, = 4.5, = 0.3, f3r = 0.02, ar = 0.12. Other parameters are 

^bins = 200, iVblock = 50, Nt = 40, lYprimal = 50000, iVdual = 4000, A^sub = 50. 


which have been developed in the last twenty or so years may be combined to provide a 
practical solution in many instances. Working entirely with Markovian problems, the key 
elements to the approach studied here are: 

• pretend that the stopping reward process Z is itself Markovian, and by discretizing Z 
onto a suitably-chosen hnite set of values we estimate the transition probabilities of 
this hnite state Markov chain by simulation ( this is the approach of Barraquand & 
Martineau 12]); 

• solve the optimal stopping problem for this hnite state Markov chain by dynamic pro¬ 
gramming; 

• use the solution to generate a stopping rule whose performance is evaluated by simula¬ 
tion; 

• use the dual characterization of the value of the problem (see [Tl|, [S], [I]) to hnd a 
hedging martingale. 

This approach is a very general methodology that delivers upper and lower bounds on the 
price that are generally reasonably close, but more importantly it provides ehective recipes 
for action. It is often said that the price a bank charges for a derivative is more to do with 
the cost of hedging that derivative than with any number that comes out of some model; and 
the approach we advocate here puts that into effect. Indeed, at each time the analysis we 
propose tells the seller of the option what hedge he should use - all he has to do is to hedge 
the approximate value at the next time step. That approximate value is a simple function of 
the stopping reward at the next time step. Similarly, the approach we use provides a compact 
solution for the buyer of the option; at each time, he calculates the value of immediate 
stopping, and stops if and only if this value is in some hnite union of intervals. 

This approach is remarkably successful in many of the examples we have studied, providing 
bounds which are often within 5% of each other. Since this sort of error is inherent in the 


18 











estimates of input parameters, or in the assumption that those parameters are constant over 
time, there is really little beneht in getting the bounds very much tighter. We would ideally 
like to have methods which (if they cannot get the derivative price exactly) can give bounds 
which are apart by, say, 1 basis point. This is an industry-standard criterion ... but where 
does it come from? Does a bank really care if their calculation of the price of a derivative 
is out by 1 cent in $100? Of course not! The Ibp criterion really comes from the desire to 
delta-hedge the option; so we want to vary the prices of the underlyings by ±1% and then hnd 
the change in the price in order to put on the delta hedge, and at this point Ibp accuracy is 
a relevant requirement. But the approach here gives us the hedging strategy by a completely 
different route - there is no delta hedging, only the hedge that comes from the dual approach! 
Moreover, getting Ibp accuracy from a simulation method is already rather over-optimistic. 

So what happens when the upper and lower bounds are further apart, as in the hxed-strike 
Bermudan Asian option? This is not problematic conceptually; the lower bound is what the 
buyer objectively thinks the option on its own is worth, the upper bound is objectively what 
the seller expects hedging this derivative on its own to cost him, but either side may move 
beyond the bound if by taking on the contract they lay off risk elsewhere in their portfolio 
- this is the most basic reason for a market in derivatives. There is nothing difficult or 
contradictory in a market where the bid price is below the ask price - this is normal! If we 
hnd a situation where the bounds are far apart, and we feel it is important to bring the 
bounds closer together, what can we do? There are four places where error enters into this 
approach: 

1. the assumption that the reward process Z is a real-valued Markov process; 

2. the error from discretizing Z values into bins, and deriving the transitions from simula¬ 
tion; 

3. simulation error in evaluating the stopping strategy; 

4. simulation error in evaluating the hedging strategy. 

We can reduce the last two by doing many more simulations, and the second can also be 
addressed by taking more bins and doing more simulations, but the first error is intrinsic - 
we can do nothing to reduce it other than change the problem in some way. The hxed-strike 
Bermudan Asian option illustrates this well, because as we discussed in Example 13.41 the state 
variable for this particular problem really has to be the two-dimensional vector {St, At), and 
the crude notion that we can get a good approximation just by looking at At on its own is 
shown by our calculations to be wide of the mark. Now of course we could work to exploit the 
specihc features of this problem to devise a problem-specihc solution (and the hnite-diherence 
calculations in Longstah & Schwartz [12] and Rogers & Shi [15] are instances), but this is 
contrary to the generic nature of the approach presented here. If we wanted to continue to 
use this approach, we might try to bin the values of the bivariate process {St, At), which is 
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of course conceptually no different from the binning of the scalar-valued reward process, but 
we may expect that the coarser binning we can expect from a two-dimensional underlying 
process will raise the second type of error. 

To summarize, then, the approach explored in this paper: 

• is completely generic - the same code gets used with only changes in the Markov process 
and the reward function; 

• always gives bounds, which are often within the range of error introduced by estimation 
or questionable modelling assumptions - and always within the proht margin required 
for an OTC product; 

• is largely insensitive to dimension; 

• gives simple and explicit exercise strategies, and hedging strategies; 

• requires only the ability to simulate a step of the underlying Markov process. 

The approach of this study puts together earlier discoveries from the last twenty years. While 
it may be premature to declare that Bermudan options in high dimensions are now a solved 
problem, what we have seen is that there is a general approach which proves to be at very 
least a good start on the problem, even if for a particular question we may want to dig deeper. 
Once we set the question in the context of estimation error and model uncertainty, it may 
indeed seem pointless to dig deeper in any case. 
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